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FOREWORD 


This is the final report of a several-years research project entitled 
"The Structure and Circulation of the Lower Martian Atmosphere," sponsored 
by the Planetary Atmospheres Section, MSA Langley Research Center, under 
Grant No. NGL-26-006-016. 

In the exploratory stages of the project, it became apparent that a 
thorough investigation of the radiative-conductive transfer of thermal 
energy in the soil and atmospheric subsystems leads to a deeper under- 
standing of the temperature field and static stability in an atmosphere 
which has a radiative responsiveness much greater than that on Earth. The 
research reported herein reflects this basic orientation. 

Data from various Mariner spacecraft were readily made available by the 
Technical Information and Documentation Division of the Jet Propulsion 
Laboratory, Pasadena, California. Dr. W. A. Baum, Director of the Planetary 
Research Center at Lowell Observatory, Flagstaff, Arizona was most helpful 
in providing photographic information on the planet Mars. Finally, the 
authors wish to express their gratitude to Drs. R. A. Hanel and J. C. Pearl 
of the Goddard Space Flight Center for generously releasing Mariner 9 
preliminary data. 
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1.0 Introduction; Objectives of Report 


Among those characteristics of the Martian atmosphere 
which have important bearing on the nature of predominant 
atmospheric processes, one of the most essential is the 
capability of the atmosphere to respond rapidly and in- 
tensely to imbalances in the atmospheric radiation field. 
Because of this capability, it is likely that the Martian 
atmospheric circulation has a significant component which 
is directly driven by thermal gradients produced by dif- 
ferential radiative heating, with the resulting motion 
field having little impact on the redistribution of the 
thermal field itself (Gierasch and Goody, 1968). 

Evidence for conditions conducive to a large radiative 
response capability on Mars has a firm basis in the data 
returned by the Mariner 4,6,7, and 9 spacecraft series. 
The Mariner 4 occultation experiments provided first in- 
dications that surface pressure was characteristically 
less than 10 mb, while spectroscopic measurements, par- 
ticularly those performed by Mariner 6 and 7, confirmed 
a high concentration of CO^ asmajor atmospheric con- 
stituent. Inversion of infrared interferometer spec- 
trometer (IRIS) measurements (Mariner 9) have provided 
representative temperature soundings, although much work 



remains to be done before the impact of suspended dust 
aerosol on IRIS-derived temperature distributions can be 
assessed . 

In view of the importance of radiative processes in 
providing the elementary drive for Martian atmospheric 
circulation of various space scales, it is essential that 
we understand in detail the radiative exchange mechanisms, 
particularly those in the near-surface atmospheric layers. 
In the planetary boundary layer, and much more so in the 
surface boundary layer, a significant diurnal variation 
in thermal structure is induced in response to large sur- 
face temperature variations. It is in these layers that a 
large fraction of the conversion of insolation into at- 
mospheric internal, potential, and kinetic energy takes 
place . 

The direction taken by the research reported here has 
been influenced greatly by the preceding considerations. 
The objectives of this report are to describe a research 
program and numerical model designed to investigate the 
role of radiative processes in controlling the subhourly 
history of the Martian atmospheric thermal structure, and 
the time-dependent interplay between physical processes 
responsible for surface heating and those determining the 
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response of the atmospheric boundary layer to that heating. 
This type of investigation can contribute to a better 
understanding of the dynamical structure of the Martian 
planetary and surface boundary layers, which is important 
for Martian lander vehicle design and mission planning, 
and to a more complete characterization of the driving 
mechanism for larger-scale circulation. 

As discussed in Section 2, the data base provided by 
the Mariner spacecraft series is sufficient to allow in- 
ferences into the basic physical mechanisms of radiative- 
conductive heat transfer in the Martian atmosphere-ground 
system. These inferences provide the basis for the numer- 
ical simulation model, which is described in Section 2.2 
and Section 3 . Section 4 relates the basic numerical ex- 
periments which have been performed with the model, and 
interpretation of results is undertaken. The theoretical 
framework for an expanded version of the model is de- 
scribed in Section 5, and Section 6 summarizes major 
findings of the research program. 
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2.0 A Physical-Mathematical Model of Martian Ground- 


Atmosphere Thermal Boundary Layers 


2.1 Basic Physical Mechanisms 

Inferences relating to the basic physical processes at 
work in the Martian ground-atmosphere system may be made 
on the basis of the observational evidence outlined in the 
Introduction. These inferences can then be used as a guide 
in developing a numerical model of heat exchange in the 
thermal boundary layer. 

Of crucial importance in this context is the evidence 
for a relatively low surface pressure (~6 mb) and high 
concentration of carbon dioxide as atmosphere constituent. 
This combination of the predominance of a radiative ly-active 
constituent gas and a relatively tenuous atmosphere results 
in an intense radiative response capability, so that heat 
exchange in the troposphere is dominated by radiative 
transfer for thermal perturbations of characteristic scale 
greater than 100 m. 

Goody and Belton (1967) have obtained the preceding result 
through calculation of the time for relaxation of thermal 
perturbations by radiation, molecular conduction and 
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turbulent convection in a dust-free carbon dioxide model 
atmosphere. It is estimated that the Martian atmosphere 
radiatively responds to surface temperature variation to 
a height of several km, particularly if the surface vari- 
ation is large over diurnal time periods. 

The intense radiative response capability has a pro- 
found impact on Martian atmospheric dynamics. Gierasch 
and Goody ( 1968 ) contend on the basis of scaling estimates 
and relaxation time calculations that to first order the 
nonlinear advective feeback of mean velocity on the tem- 
perature field may be neglected in the Martian atmosphere. 
Hence, it is likely that a model which neglects advective 
temperature changes due to the mean wind field can provide 
meaningful information on the time-variation of thermal 
structure in the lower atmosphere (Gierasch, 1971 ). 

A second crucial body of information provided by the 
Mariner observations is the evidence for a global extent 
of sand-like Martian soil material and lack of oceanic 
areas. These features permit large ly-unattenuated in- 
solation to induce a large diurnal soil surface temperature 
variation, especially in the absence of suspended dust 
particulates. Both Earth-based and Mariner radiometric 
scans indicate a diurnal variation of brightness temperature 
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of 50-80K at midlatitude s . In view of its radiative 
responsiveness, we may expect that the Martian atmo- 
sphere will respond significantly to this surface tem- 
perature variation, in addition to its direct response 
to solar heating. The latter is especially important 
if there exists an appreciable amount of suspended dust 
particulates. 

Details of the coupling between Martian soil and 
atmospheric temperature variation due to radiative ex- 
change depend crucially on the structure of the C0 2 ab- 
sorption bands. Since most of the 15/* band absorbs strongly, 
even under the reduced surface pressure conditions of Mars 
(~6 mb), a portion of the coupling will occur through the 
near-surface and interface layers, where radiative transfer 
will approximate a diffusion process. On the other hand, 
spectral "gappiness" and weak-line transmission, each 
enhanced by atmospheric tenuity, will permit some long- 
range exchange of energy from the planetary surface to 
several scale heights and towards space. Absorption of 
direct solar radiation in the 4.3/* band, and the presence of 
lines of intermediate intensity, further complicate the 
heat transfer by radiative exchange. 
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The preceding considerations point to the necessity of 
considering the thermal coupling between ground and atmo- 
sphere subsystems in any model of the diurnal variation of 
Martian thermal boundary layer structure. Molecular thermal 
conduction in the soil material can account for an ap- 
preciable fraction of the net heat flux balance at the 
soil-atmosphere interface (Pallmann and Dannevik, 1972; 
Sellers, 1965), and requires knowledge of the subsurface 
temperature profile for its determination. Hence, cal- 
culation of the temporal variation of boundary layer ther- 
mal structure requires a simultaneous forward time-inte- 
gration of the subsurface temperature distribution. 
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2.2 Development of the Mathematical Model 


2.21 Physical Setting 

The atmospheric model is taken to be that of a pure 
carbon dioxide atmosphere of 50 km depth, in local therm- 
odynamic equilibrium, and exhibiting plane-parallel strati- 
fication. (However, as discussed later, the impact of sus- 
pended dust particulates on gaseous absorption will be taken 
into account by the model). Doppler broadening of spectral 
lines is assumed negligible, as well as the effects of 
molecular scattering. 

The lower solid surface is taken to be either a flat, 
uniform sand-like material, or solid carbon dioxide.- The 
surface diffusely reflects solar radiation in accordance 
with a radiometric albedo, and emits as a grey body. In 
the case of solid carbon dioxide, the interface temperature 
is assumed fixed at the solid-vapor equilibrium sublimation 
point, and the subsurface temperature field is isothermal. 
Por all other cases, the ground-atmosphere interface 
temperature, as well as the subsurface profile, is deter- 
mined as a function of time by the model. 
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2.22 Governing Relations 


For the preceding physical setting, the applicable 
heating-rate' equation takes the following form: 

, o<xo, s ,i>o 

I U) 

In (1), x is the vertical coordinate, measured from the 
top of the model atmosphere (x=0), x g is the ground- 
atmosphere interface level, and x d is the lowest subsurface 
level under consideration. Other notation is defined in the 
Glossary of Symbols. 

The preceding equation neglects the production of heat 
by viscous dissipation, the aavection of temperature by 
a mean velocity field, and the contribution by the rate 
of working through shear stresses. These approximations 
are complementary to the assumption that the advective 
feedback of velocity on the temperature field is neglible. 

In addition, it is assumed that no phase changes ocurr 
in the regions away from the interface level. 

The right-hand side of (1) represents the convergence 
of total heat flux at the level x and time-t . The first 
term represents the contribution of molecular thermal 
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conduction, and the second the impact of radiative flux 
at all wavelengths. Formulation of the preceding quantity 
is the subject of an earlier report (Pallmann, 1968). 

The spectral radiative net flux for the case of local 

thermodynamic equilibrium may be expressed as 

o x , . * s . 

5 feM ) e <^<4 {2) 

0 i 

In this expression, S v (0) is the intensity of radiation 
incident at the top of the atmosphere, and G y (x s ) is the 
upwardly-directed radiation at the soil-atmosphere inter- 
face level. The planetary radiation field is assumed to be 
axially symmetric. It may be seen from (2) that the total 
spectral radiative flux is a functional of the atmospheric 
temperature distribution, through its dependence on the 
Planck function. 

To complete the mathematical model, boundary and initial 
conditions are required for (1). 
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2.23 Boundary and Initial Conditions 

At the interface between atmospheric and soil subsystems, 
a discontinuity in physieo-opt ical properties occurs, so 
that (1) is not directly applicable at this level. We may, 
however, derive an internal boundary condition on the as- 
sumption that the total heat flux must remain continuous 
across the interface. This is equivilent to the requirement 
that there can be no net accumulation of thermal energy in 
a layer of vanishing thickness, since this would lead to 
an unbounded temperature change. Considering the fluxes 
of solar and planetary radiation and molecular-conductive 
heat flux in the soil and atmospheric layers adjacent to 
the interface, the requirement of continuity of heat flux 

at x=x takes the form 
s 

X O CO % 

0 J | J o 

o to j 

-»0 v 1 0 

In this relation, the heat flux due to phase changes 
has- been omitted and solar radiation incident at the sur- 
face is assumed diffusely reflected. 

To (3) may be added the requirement that soil and atmo- 
spheric temperatures are equal at the interface level. 
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This constitutes an assumption that temperature is con- 
tinuous across the interface, which is complementary to 
the requirement of local thermodynamic equilibrium for 
the atmospheric subsystem. Thus, we have 

(i)) 

in addition to relation (3). 

In the case that the lower surface is solid carbon 
dioxide, the requirements (3) and ( ^ ) may be simplified 
since the interface temperature may be assumed equal to 
the CO^ solid-vapor equilibrium value. Hence 

* TlX*-®,*) = T sub 

where T Sub is the sublimation temperature for carbon 
dioxide at the surface pressure. 

At the upper and lower bounds of the solution domain, 
the temperature is assumed fixed at a known value: 

T(p,«=T 0 .T.Wj^zT, 

(o) 

If the absorption coefficients, thermal conductivity, 
density distribution, and physical properties at the 
interface are assumed given, then a known initial temper- 
ature distribution is all that is required to integrate the 
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system (1), (3) and (4) or (5), and (6) forward in time. 

The complete system of equations is presented in the fol- 
lowing subsection, followed by some analytical simplifications 
to increase its tractability . 


2.24 Formulation of the Mathematical Problem and 
Analytical Simplications 


The relations presented in the preceding discussion 
may be combined to pose the following boundary-initial- 
value problem: 


Governing equation: 

eC ft = + (r3 1 f C S j ^ ^ s^’ f 

' 3^,t) cy pM fi ile^l-U^cAsA’) <?K V ( ^ ^ 

' P X* “ 



X 5 < * < > O 


(7b) 


Boundary conditions: 

At x=x : 
s 

k s -£^T 5 4 (X s -^» t -)+2.uS i $06,i)ex 

O ct> - 

+ 2 * 7? ^ f \ eJCdfe/U 1 ~ ° , y 

or, 

T s (X 5 ^t) = T U 5 ^,t.^T SuV ,. 
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At x=0: 


T(»,t1 = T 0 , t»° 


At x=x . : 
a 


(9a) 


V*A,t\*T A 


i^/O 


(9b) 


Initial conditions: 

T s.tO 

TU.ol = T L tO 

The system (7) through (10) is sufficient to determine a 
unique solution for the distribution of temperature through- 
out the region for any time -t>o . 

Before adapting the system (7) througn (10) for numerical 
analysis, some simplications in the analytical formulation 
can be made without seriously affecting the solution. 

First, ;we consider the possibility of performing im- 
plicitly the integration over the cosine of zenith angle 
indicated in eqs. (7) througn (8). Consider, for example, 
the integral in (7) having the form 

14 


, 0<x<x 5 (10) 



( 11 ) 


• i, i 

} L~\ 

This expression involves the exponential integral of 
order n, defined as (Sparrow and Cess, 1966) 

E„im (12 > 

To a sufficient degree of approximation, the integral 
of order 3 behaves as a simple negative exponential: 


a Pin* e b) \ b=i.<,4 (13) 

3 J 

This approximation has been discussed by Elsasser (I960), 
Goody (1964), and Kondratyev (1969). Using the formula 


dr 


"= z >, 


• 4 • 


the integrals involving E„ may be replaced by approximations 
similar to ( 13 ), and the explicit integration over zenith 
angle is eliminated. The expression (2) for spectral net 
flux then appears as 
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F v (x,-0 = zirS^o^ +i7rb' 1 G i/ 0f5)^p(->o^x y(< is) 

-*•*' fcll.O A[e, f 1-bj‘t^J.')] Jf 

-Zub"' ( S^.O^pW J? (14) 

where ^ now refers to the cosine of solar zenith angle at 
the time t. 

Similarly, the internal boundary condition (8) may be 
rewritten as 

ks E i - k ^' £rT *V^WaFr6; S ,-t)e,pt-ibeK,cis')j v d s 

Op * 

CO *5 j / I r \ 

+zn ^0-tf 5 )S le> )«*p(-^ ) j a ^ v ~ ° 

p v p ' 

It will be convenient for the computational version of (15) 
to simplify the grey body emission term frT^ . This term 
may be expanded around a reference temperature T s<j , which 
will later be identified with the surface temperature at 
the beginning of a timestep in the numerical integration. 
Thus, w e may write 

£ <fT/w 5 ~0,-t)- H6(TT S * [T 5 W5-^)-^Ts 0 l (16) 

Since, on a relative scale, the surface temperature change 
over a single timestep is small, the representation ( 1 6 ) is 
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sufficiently accurate. With the approximation (16), 
the finite-difference form of ( 15 ) appears as a linear 
algebraic equation for the temperature at the 

end of a timestep. 


3.0 Design of the Numerical Algorithm 

3.1 Integration over Frequency and Computational 
Form of the Transmission Function. 

The integration over frequency for the radiative flux 
terms in (7) and (8) will be accomplished by dividing the 
CO,-, absorption spectrum into small increments, within which 
the Planck function and solar spectral intensity will be 
assumed constant. This procedure leads to the consideration 
of the CO 2 transmission function, defined by 

vj+tr.- s 

t V - Av ' 

*• i 

where is the frequency interval of interest, and s^ 

and S 2 are the endpoints of the geometric slant path. It 
should be noted that the mass absorption coefficient 
is a function of temperature and equivalent pressure along 
the path length, and the transmission function must take 
into account the variation of these gas-kinetic properties 
over the distance IS, - S 1 1 

The intricate fine structure of the C0 2 vibration- 
rotation bands make the CC >2 absorption spectrum suitable 
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for modeling with the quasi-random transmission function 
(Goody, 1964). The essence of this model is that lines 
are placed at random with respect to frequency, implying a 
Poisson distribution of the spacing between neighboring 
lines. The well-developed theory of single-line models 
is then readily utilized to construct a band model with 
overlapping lines. The analytical form of the resulting 
model is given by 

Ss.SflVSi/UTrap*) 

{b = 2D-tf/d (17) 

where: 5 . -.average line intensity in 

the i-th frequency interval 

« mean half width 

c):mean line spacing 

Tpl, '.Bessel functions of the 

first kind of imaginary 
argument . 

Prabhakara and Hogan (1965) have utilized the laboratory 
data of Stull, Wyatt and Plass (1963) to construct a re- 
presentation of ( 17 ) for the carbon dioxide spectrum, 
assuming that gas-kinetic properties do not vary along 
the geometric path. The tabulated band parameters are. 
given by c,= S 0 p 0 /(ZTT<* 0 ) , c 2 = nr« 0 /( ?0< i^ . 


19 



With these parameters, the arguments of the- 
function (17) may be calculated for non-standard pres 
sure and temperature through the relations 

e.pjf 

The tabulated coefficient P accounts for the dependence 
of line intensity on temperature in the 15/i-band. 

From ( 17 ) we way derive the following limiting forms: 

so— *• «*pl- f»^) 

v (18) 

c *p Z ]j S >3 

These are the weak- and strong-line limits, respectively. 

It must be stressed that (17) is valid for a path length 
of uniform gas-kinetic conditions. In the governing equa- 
tions (7) and boundary conditions (8), the transmission 
function must be evaluated over extended paith lengths, so 
that variations in pressure and temperature are signif- 
icant. In the numerical model, this effect is approx- 
imated by dividing a given extended path into a number of 
shorter segments, within whicn the gas-kinetic conditions 
do not vary appreciably. The transmission functions for 
each of these segments are then calculated, and their 
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product is used to approximate the transmission of the 
extended path. 

For a given spectral interval and geometric path 
length, either the weak- or strong-line forms of (17) is 
used, depending on the value of & for the particular cal- 
culation. 

In those atmospheric strata assumed to be dust-laden, 
the gaseous transmission is reduced by a factor expC-p^ax), 
where (3^ is an "effective absorption" coefficient due to 
dust, and Ax is the thickness of the layer. 

3.2 Finite-Dif ference Form of the Governing System 
of Equations 

The system (7) through (10) must be put into a form 
suitable for computational analysis. The basic procedure 
is to approximate the integrals in (7) and (8) by quad- 
rature formulae, and the differentials by finite dif- 
ferences. It should be recalled that' integration over 
the zenith angle has been performed implicitly by in- 
troduction of the "aiffusivity factor" b (eq. ( 13 )). 
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As discussed in Section 3*1* integration over frequency 
is approximated by summation over a number of frequency 
intervals embracing the 1-6 and 12-18^ vibration-rotation 
bands of carbon dioxide. The CC^ band parameters and 
raw solar spectral intensities are presented by Pallmann 
(1968). 62 increments of approximate width 50 cm“^ are 

utilized, and the lower surface is assumed to be a "grey" 
emitter. 


With the i-subscript referring to a frequency increment, 
(7a) and (8) take the form 

l (K^ + iT\l[Sjo)Tu>i , )-fG,.u 5 )TlUx-x 5 )] 

2>t <?x l i u *• ' u u 


k*g‘- kg ) + zn?[^ 5 e' 6.i<V)^[b(s-y s >] 4s' 

+ ci-ois) S (x*/0 * 0 


(19) 


( 20 ) 


Each of the transmission functions is approximated by a 
product, as indicated previously. For example, we may 
write 


< T.u,-^) = Tr ) 

i=. L 1 •- 


( 21 ) 


wnere the product index j extends over all J intervening 

segments (x. -x. ) between the endpoints of depth x, ,x . 

\ i— 1 ± 2 
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Integration over the depth coordinate in (19) and (20) 
is represented by a summation over a number of atmospheric 
strata. For convenience, these strata are taken as the 
same ones used in the product algorithm (21). Using the 
subscript index m for a given level between the surface 
and 50 km, (19) and (20). assume the form 


M . M ( 22 ) 

’** w -r* 'iyr, J 


(23) 


Finally, the derivatives in (22) and (23) are replaced 
by finite differences. A simple forward difference is used 
for the time derivative, while modified centered differ- 
ences are employed in the space derivatives. The re- , 
suiting equations, including the product representation (21), 
are quite lengthy and will not be reproduced here. The 
corresponding finite-difference form of (7b) is 


h+i _ 1 

I*-- I 




e~ c = * -T+i" } (2 4) 

— ^ (AX) 1 m ‘" S" 1 S* + J ) A v ' 
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where the superscript refers to the timestep, and m 

d 

is the lowest soil level under consideration. In (2*0, as 
well as the finite-difference form of (22) and ( 23 ), the 
Duf ort-Frankel two-level explicit differencing scheme 
(Ricntmyer and Morton, 1967) has been employed. This 
scheme is unconditionally computationally stable. 


From (15) and (16), the finite-difference form of the 


internal boundary condition takes the form 

i T >1+l T"* 1 m*' m+i _3k . m. •, * N 

k. laide. -kkln -Hfrlj (T 5 - *T S ) + 

A <s " " * 

l IAA r i I ^ 


l m - 


Note that the soil surface temperature for the previous 
Cn-ih’) timestep is used for the expansion (16). Again, 
the product algorithm (21) is used to evaluate each trans 
mission function in (25). Each of the factors in the 
product is evaluated using the temperature distribution 
from the previous timestep. 

We now turn to the question of the size of space and 
time increments to be used in the simulation. In view 
of the weak-and strong-line radiative transfer, we can 
expect a relatively complex atmospheric temperature pro- 
file, especially near to the surface, where large 
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temperature gradients will result in significant molecular 
thermal conduction. A reasonably high vertical resolution 
in the finite differencing is required to accomodate this 
feature, along with the requirement that temperature and 
pressure do not vary substantially over a given layer, so 
that the product algorithm (21) for the transmission func- 
tion remains accurate. We have chosen 7 layers approximate- 
ly logarithmically spaced in the lowest 1 km, and generally 
a 1 km spacing in the region between 1 and 50 km. A 2-cm 
increment was used for the subsoil layers. 

From the estimates of Goody and Belton ( 1967 ), the 
characteristic time for relaxation of thermal perturbations 
having scales of the order of several meters is in the 
range of 10 J sec for the Martian troposphere. On these 
physical grounds, we may infer that a computational time- 
step of similar size is required to adaquately resolve in 
time the propagation of temperature waves in the near- 
surface layers. A 15-min time-step was used for all 
simulations reported in Section 4 . 
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3.3 Input Data and Computational Procedure 

The parameters and input data required for a given 
simulation are as follows: 

(1) initial distribution of atmospheric temperature, 
density, and subsoil temperature; 

(2) soil conductivity, density, and heat capacity, 
each assumed constant; 

(3) atmospheric conductivity and heat capacity at 
constant pressure; 

(4) soil surface albedo and emissivity; 

(5) C0 2 band parameters and solar spectral radi- 
ances for each frequency interval; and 

(6) orbital parameters, relating to the Martian day 
under consideration, and latitude of simulation. 

The time-variation of solar radiation entering the top 
of the atmosphere is computed from the relation presented, 
e.g., by Sellers ( 1 9 6 ) . 

Thermal parameters of the soil material were derived 
through a parametric procedure, in which soil thermal in- 
ertia was adjusted until the amplitude of diurnal surface 
temperature variation matched that suggested by Mariner 
radiometric scans. Numerical values of input parameters 
for each simulation reported will be presented in Section 4. 
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The computational procedure followed for each simula- 
tion begins with the calculation of transmission functions, 
solar spectral radiance, and atmospheric spectral emission 
corresponding to the initial time and temperature dis- 
tribution. These quantities are then combined to form the 
total (conductive-radiative) heating/cooling function 
represented by the right-hand-side of (19). 

In addition, the downward radiative flux incident at 
the surface is calculated, as required in (25). Eqs. (9)> 
(24) and (25) are then solved for the subsoil and surface 
temperature at the end of the first timestep, and the new 
distribution of atmospheric temperature is computed from 
(23), with the boundary condition (9a). This procedure 
generates the subsoil and atmospheric temperature distri- 
bution for the end of the first timestep. New solar radi- 
ances are then calculated for the advanced local time, and 
the new temperature field is used as initial data for the 
execution of the next timestep. In this manner, the pro- 
cedure continues through the total interval of time desired 
for the simulation. In most numerical experiments, this 
interval corresponds to 2k Martian hours, or 96 computa- 
tional timesteps. 
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In the cases where the lower solid surface is assumed 
to be solid carbon dioxide, the internal boundary con- ■ 
dition (25) is replaced by the requirement that temper- 
ature remains fixed at the CO^ - sublimation point. The 
subsoil temperature field is assumed isothermal, so (2k) 
is no longer required. 

In all, there are 82 computational levels and 62 fre- 
quency intervals considered for each timestep. On the 
CDC 3300, approximately 20 sec of computing time are 
required for each step in the forward time-integration. 
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4.0 Basic Numerical Experiments and Interpretation 


of Results 


4.1 Description of Experiments 

The results contained in this chapter are based on a 
series of time-dependent numerical simulations of the 
Martian thermal energetics. Emphasis is given to the 
effect of a uniform dust layer on the transfer of radi- 
ation. This scientific interest evolved from the ob- 
servational evidence of Mariner 9* ?or the simulations 
considered, a comparison is made between the results for 
a pure C0 2 atmosphere and those for a partially dust-laden 
atmosphere. The midi at itudinal and polar conditions are 
parametrically represented in table I. 

TABLE I 


Parameter Mldlatitudinal-Experiments S. Polar Region Experiments 

(Mariner 7) (Mariner 9) (Mariner 9) 


Initial data 

S-'oand Mar 7 

IRIS 

IRIS 

IRIS 

Initial time 

0300 hr 

1900 hr 

1300 hr 

1800 hr 

Season 

N.H. Autumn 

S . H . Sum . 

S.H. Sura. 

S.il. Sum 

Latitude 

+ 33° 

-38° 

-80° 

-80° 

Declination 

-6° 

-2 3° 

-23° 

-23° 

Surface aloeao 
Surface emls- 

0.25 

0.25 

0.95 

0.30 

sivlty 

0.90 

1.0 

0.90 

1.0 

Thermal con- 





auctivity, 
CO, gas 
(css) 

1.33 x 10 3 

1.33 x 10 3 

1.33 x 10 3 

1.33 x 10 3 

Soli thermal 


C 



inertia 

(cgs) 

5.90 x 10 3 

5.90 x 10 3 


5.90 x 10 3 
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The midlatitude simulations involve a dual structure 
in which the ground-atmosphere subsystems are coupled 
by an interface energy balance, thereby accounting for 
temperature variations in the soil. For the polar ice 
cap simulations, the temperature of the ice cap surface 
is held fixed at the sublimation temperature of CC^. 

An additional simulation is reported for a frost-free 
lower boundary. 


4.2 Results and Interpretation 

In this section, the material will be divided into 
three categories, based on the initial input data used 
in the simulations. Due to the large volume of reliable 
information derived from Mariner 9, most of the statements 
to follow pertain to simulation results based primarily 
on Mariner 9 data. 


4.21 Mariner 7 Midlati tudinal Experiment 

Numerical results of the first simulation evolved from 
Mariner 7 initial input data, i.e., temperature & pressure 
as a function of altitude derived from the S-bana occul- 
tation experiment. Calculated atmospheric temperature 
distributions in the lowest 400 m, at specified times during 
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the Martian day, are portrayed in Pig. 1. The curvature of 
these profiles indicates that molecular conduction produces 
cooling at 0400 LMT , warming at 0900, no change at 1400, 
and cooling once again after 1800 LMT. (Note the logar- 
ithmic height scale). A dual structure 



Fig. 1 Calculated temperature 
profiles in the lowest 
400 m for various hours 
of the Martian Day. 
Data: Mariner 7, mid- 

latitude (see Table I). 


in the diurnal temperature wave propagation .pattern is 
evident. One wave, primarily confined to the lowest 150 m, 
is superimposed on a second wave of larger characteristic 
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scale. Goody and Belton (1967) have predicted such a 
structure, based on the difference between tine of re- 
laxation of thermal perturbations by radiative and con- 
ductive processes in the Martian atmosphere. Their 
analysis indicates characteristic scales of 15 m and 1.5 km, 
respectively, for the molecular-conductive and radiative 
wa ve s . 



Fig. 2 Time-variation of molecular- 

conductive heating rate at the 
12.5 in computational level, cor- 
responding to the lowest 25 m 
layer. 


Fig. 2 reveals the relative significance of atmospheric 
molecular conduction. At 1800 LMT, the conductive heating 
rate at the 12.5 m computational level reached -3.9 x 10 J K 
sec” 1 (-3.4K day” 1 ), representing 22% of the total 
(radiative-conductive) heating/cooling function. 
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Fig. 3 Atmospheric and soil temperature profiles below 
300 m (3 hm), showing evidence of dual-wave 
structure in thermal boundary layer. The subsoil 
temperature below 60 cm of depth is one of sev- 
eral simulation temperatures chosen. 

In Fig. 3, the typical diffusion-like upward propagation 
of the perturbation at 0900 LMT may be traced through the 
profile for 1^00 LMT when a second disturbance has been 
initiated by surface cooling. Again, an additional wave of 
larger scale is indicated, extending well into the region 
above 300 m. 


H.22 Mariner 9 Midlatitudinal Experiment 


The second simulation was conducted with the purpose of 
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investigating the effect of dust an the Martian thermal 
energetics. The results presented are based on two 
unique 24 hour simulative investigations. 

The first case, designated "dust-free" utilizes, as 
initial data, the Mariner 9 IRIS inversion profile men- 
tioned previously, with a local starting time of 1900 
LMT. The second case used an identical initial sounding 
and starting time, but includes a simulated ground-based 
dust layer of 30 km thickness; this case will be termed 
"dust-laden" . 

Since over the night-time hours the radiative-con- 
ductive model is thermally dissipative, minor inconsisten- 
cies between the arbitrary subsoil and surface temperature 
and the atmospheric sounding tend to disappear over the night- 
time hours so that by sunrise, the entire temperature profile 
is thermally "initialized". 

Fig. 4 depicts modeled atmospheric soundings for 0500 
and 1600 LMT for each simulation. With regard to the 
dust-free case, basic regimes are evident. The lowest 
regime, from the surface to about 4 km, includes those 
atmospheric strata dominated by strong-line, short-path, 
diffusion-like radiative exchange associated with soil 
surface temperature variation. (This region actually 




(a) (b) 

Fig. 4 (a) Modeled vertical temperature profiles at 

0500 and 1600 local Martian time (LMT), for 
the dust-free (solid line) and dust-laden 
(dashed line) conditions. Also shovm is the 
initial (1900 LMT) profile obtained from in- 
version (Hanel, et al., 1972) of Mariner 9 
IRIS measurements by assuming a dust-free 
atmosphere. (b.) profiles of the difference 
AT = T - T, in temperature between the 
dust-free ( subscript : c ) and dust-laden (sub- 
script :d) simulations, at 0500 and 1600 LMT. 
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exhibits a further fine structure, to be elucidated in con- 
junction with Fig. 2, which may be attributable to the 
diffusion-like dual action of molecular conduction and 
strong-line radiative transfer). The region between 4 and 
31 km is directly affected by surface temperature variation 
primarily through weak-line exchange, but also undergoes 
direct solar heating and weak-line losses to outer space. 

The regime above 31 km is dominated by the upper boundary 
conditions, i.e., direct solar heating and emission to 
outer space, and is affected little by soil surface temper- 
ature variation, which lags solar forcing by roughly 4 
hours . 

Soundings of the dust-laden case exhibit essential 
differences from those of the dust-free simulation primarily 
in the middle regime (4-31 km). Over the daylight hours, 
this region approaches an isothermal stratification, in 
agreement with recent Mariner 9 findings showing a tendency 
towards this structure in the dust-laden mission phases. 

The top of the dust layer acts as an "effective radiation 
surface", similar in its radiative impact to the solid 
soil surface. Diurnal temperature variation at the top of 
the dust layer is nearly 50 % of the amplitude of surface 
temperature variation. Thus, the upper regime (above 31 
km) responds to this effective radiation table in the same 
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way as the lowest regime responds to surface temperature 
variation in the dust-free case. A statically-unstable 
layer is generated at the top of the dust layer in the 
afternoon hours. 

The initial (1900 LMT) Mariner 9 IRIS profile was derived 
by inverting the IRIS radiance measurements obtained during 
a dust-laden mission phase, to generate an emission profile, 
assuming a dust-free atmosphere. As indicated in Pig. 4, the 
profile approximates a smoothed version of the modeled 1600 
LMT output for the dust-laden simulation, with decreasing 
fidelity towards the surface. 

The curves in the right-hand portion of Fig. 4 show pro- 
files of the difference in temperature between the dust- 
laden and dust-free simulations, at 0500 and 1600 LMT. 

At 0500, the dust-laden profile is warmer than the dust-free 
in the lowest 400 m, an indication of the shielding effect 
of the dust blanket keeping surface temperature slightly 
higher in the nocturnal hours. However, in the bulk of 
tne dust layer, a cooling has occurred by 0500 LMT which 
is substantially larger than in the dust-free case. This 
is attributable to the simulated "grey-body" emission of 
the dust particles, which results in more non-absorbed 
photons being emitted per unit volume from a dust-laden 
layer than from the same layer under dust-free conditions. 
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Atmospheric and subsurface temperature profiles in the 
layers near the interface are shown in Fig. 5, for the dust- 
free case. The basic structure of the subsurface profiles 
reflects the relatively large soil thermal conductivity 
(2.3 x lO^ erg cm"l sec”-'- K”-'-) used in the simulation, 
modeling low-porosity silica. 



Fig. 5 Vertical temperature distributions in the sub- 
soil and lower atmospheric layers for selected 
hours of the Martian day, from the dust-free 
simulation. For environmental parameters used, 
see table I. 
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The atmospheric profiles exhibit diffusion-like heat 
transport characteristics, due to the dual action of 
strong-line radiative and molecular-conductive transfer. 
An "effective radiative-conductive diffusivity" may be 
estimated by li'/'t v/here L and 't are the characteristic 
length and time -scales of the diffusion process. This 
diffusivity is found to be ^ 2 x 10^ cm^ sec - -*-, on 
the basis of data represented in Pig. 5. This compares 
well with Gierasch and Goody’s (1968) scaling estimate of 
105 cm2 sec -1 , and is similar in order-of-magnitude to a 
free-convective eddy diffusivity under mildly-unstable 
conditions. This also tends to corroborate Goody and 
Belton’s (1967) estimate that the time for the relaxation 
of thermal perturbations by radiation and mild turbulent 
convection should be comparable for a perturbation length 
scale of hundreds of meters. 

In Fig. 6, each of the various heat flux components 
acting at the soil-atmosphere interface are shown as a 
function of the time of day. Short-term oscillations 
of decreasing amplitude, evident in the first two hours 
of the simulation, represent the dissipative smoothing 
of inconsistencies in the arbitrary initial temperature 
distribution in the interface and subsurface layers. 
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Fig. 6 Diurnal history of the various com- 
ponent heat fluxes at the soil 
atmosphere interface. EMIS . repre sent s 
the flux of long-wave radiation emit- 
ted by the surface, SCON the molecular- 
conductive heat flux in the soil mate- 
rial, SOL the absorbed insolation, ACR 
the absorbed flux of atmospheric counter- 
radiation, and NET the algebraic sum of 
all component heat fluxes. Solid lines 
depict results for dust-free conditions, 
while the dashed lines (shown for Ei-IIS, 
SOL, and ACR only) relate to the sim- 
ulated dust-laden conditions. Influence 
of C 02 ~gas thermal conduction is of 
lesser significance at the interface, 
and curve is not shown. 


The fact that the model does not reproduce after a 
cycle of 24 hours the initial flux values can be traced 
to at least two possible causes. The most obvious is the 
lack of the model to account for removal of heat from the 
interface layers by mechanical and/or thermal convection 



after a statically-unstable temperature distribution has 
been generated by surface radiative heating. However, it 
is also possible that the initial atmospheric temperature 
profile may not be representative of the local time and 
season relevant to the IRIS measurement on which the 
inverted sounding is based. The preliminary Mariner 9 
sounding utilized was based on an inversion method which 
assumes a dust-free atmosphere. 

Molecular-conductive heat flux in the soil material is 
a significant fraction of the net flux available for 
temperature changes, at all times during the 2k hr sim- 
ulated period. The conductive flux reaches a maximum 
about 1.5 hours before local noon, in agreement with mea- 
surements under comparable terrestrial physical settings 
(Sellers, 1965). 

Planetary atmospheric "counter-radiation" varies little 
over the diurnal cycle. Since the concentration distribution 
of radiatively-active gases remains constant in the model, 
any temporal variation in counter-radiation can be traced 
directly to variations in the temperature distribution, 
for the dust-free case. The addition of a simulated dust 
layer, assumed a grey emitter in the CO ^ absorption bands, 
results in an increased counter-radiation for all times 



of the day, as evidenced in. Fig. 6. This phenomenon accounts 
for a higher minimum morning temperature of . 198 K for the 
dust-laden case, as. compared with 193K for the dust-free 
simulation. 

Height-time cross sections of the temperature field 
in the lowest 18 km are shown in Fig. 7 . Qualitatively, 
the pattern is nearly unchanged by the presence of dust, 
in the lower 4 km. As discussed earlier, this region 
responds primarily to changes in the. interface temper- 
ature, and variation of the latter quantity, is not 
significantly modified by the presence of the dust layer 
(cf. Fig. 6 , EMIS-curves). 

The region above 4 km is modified drastically by the 
dust layer during the daylight hours. This regime 
responds to direct solar heating and, to a lesser extent, 
to interface-temperature variation. The dust-laden layers 
above 7 km develop a nearly isothermal stratification 
in the afternoon hours, in agreement with recent . Mariner 9 
observations (JPL, 1972). 


H2 



Altitude, km Altitude, km 




Local Martian Time, hr 

<t>) 

Fig. 7 . Height-time cross-sections of the modeled 
temperature field (deg. K) in the lowest 
18 km for (a) dust-free, and (b) dust-laden 
conditions. Approximate extremal lines are 
dashed. 
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4.23 Mariner 9 South-Polar Region Experiment 

The final portion of Section 4 consists of the 
simulative results for the Martian polar region. The 
essential characteristics of the atmosphere and ground 
environment are displayed in Fig. 8. Due to the tenuit 
of the Martian atmosphere, incoming solar radiation may 
be considered as beam radiation diffusely reflected by 


DUST-FREE C0 2 
ATMOSPHERE 



SOLAR 

RADIATION 


SURFACE TEMPs 148 *K 
(FROST POINT TEMP.) i / SURFACE 

A ALBEDO: 6.45 

y/ vr/ ri rrrn i / / n / n // 

ISOTHERMAL 

SUBSTRUCTURE 



Fig. 8 Schematic physical setting of the south 
polar ice cap region. Hatched area de- 
picts extent of dust layer from surface 
to 8 km. Solar radiation is diffusely 
reflected at the ground; albedo value 
of .45 represents condition of C0 ? frost 
surface in interaction with solar spectrum 
up to 5yU • 



the surface. At the dry ice surface, the frost point temper- 
ature of CO^ is the equilibrium temperature for the solid 
and gaseous phases. Therefore, the thermal condition at the 
interface was permanently fixed at 148K with an isothermal 
substructure. Further details are presented in Pallmann 
and Frisella (1972). 

Mariner 9 measurements support the assumption that a dust 
layer of moderate optical thickness extended from the ground 
to about 8 km. In the model a value of the effective in- 
frared optical thickness due to dust =1.6 was as- 

sumed for dust-laden simulations. The initial input data 
was derived from the inverted Mariner 9 IRIS data (temper- 
ature and pressure), pertinent to 1300 LMT over the south 
polar cap region. This temperature distribution together 
with the thermal structures 24 hours later, for both clear 
and dusty cases, are portrayed in Fig. 9. These distributions 
have been plotted only up to 15 km. Several noteworthy 
characteristics can be commented upon, foremost is the over- 
all cooling at all levels. In the initial sounding, 12.5 km 
corresponds to the height of the maximum temperature occur- 
rence; however, 24 hours later for both the dust and clear 
distributions the maximum temperature of the sounding occurred 
at the 7*5 km level. Below 8 km, the effect of the presence 
of dust seems evident. At first glance, a maximum cooling of 
approximately 11K in the region between 1.5 and 3*5 km is 
discernible. 
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Fig. 9 Temperature profiles for Day 1 (0-24 hr) simu- 
lation. Solid line depicts the initial Mariner 9 
IRIS inverted data, the Mash-dot’ represents the 
prognostic 24 hours calculated profile for a pure 
CO„ atmosphere, and the 'dash -dot-dot’ portrays the 
corresponding profile for a dust-laden simulation. 
All curves converge to 148K at the surface. 

Upon closer examination of the output, it became clear 
that the cumulative source layer radiative flux was larger 
in magnitude than the incoming solar flux, as can be seen in 
Fig. 10. This indicated that the temperature profiles were 
adjusting themselves toward a quasi-oscillatory thermally- 
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Fig. 10 Energy flux as a function of time (LMT) at 

the top of the dust layer (level 40). Source 
layer up, source layer down, and direct solar 
represent the upwardly- and downwardly-directed 
cumulative atmospheric radiative fluxes and 
downward direct solar radiative flux, res- 
pectively, at level 40. 

balanced distribution. The simulation was therefore extended 
out in time another 24 hours, utilizing the final sounding 
of Day 1, as initial input data.' Fig. 11 represents the re- 
sulting temperature distributions for the second day of sim- 
ulation. Once again, the maximum temperature occurred at 
7.5 km. Most importantly, a thermally balanced profile was 
established, as evidenced by the warming between the 4 and 
8 km, comparing the 0600 LMT and the final 1300 LMT outputs. 
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Fig. 11 Temperature profiles of Day 2 (24-48 hr) for the 
dust-laden simulation. Solid line depicts the 
final temperature sounding from Day 1 simulation 
used as initial data t.; also depicted are the 
1800, 2400, 0600 LMT aAd the final (t ± + 2400) 
soundings . 

The simulation was extended another 24 hours as depicted 
in Fig. 12. Comparing the initial IRIS inverted, sounding 
with the 72 hour simulated temperature distribution reveals 
a substantial cooling of about 35K at all level's. Also indi- 
cated is the existence of a super-adiabatic layer between 
7.5 and 8.5 km. It should be noted that the integral in- 
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Fig. 12 Temperature profiles at Day 3 (48-72 hr), for 
the dust-laden simulation, which started at 
1300 LMT. 

version of the Mariner 9 IRIS signals was based on an _algo- 
ritnm for dust-free conditions. This could explain part of 
the simulated cooling. It is also possible that another 
non-simulated physical transfer mechanism could reduce the 
degree of cooling. 
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The constant ice cap temperature constraint imposed in 
our simulation permits the existence of a radiative sink 
phenomenon. The emission by the warmer atmospheric layers, 
varying between' 1.7 x 10^ and 1.3 x 10^ erg cm” 2 sec 
for Day 1, initiates a downward heat flux onto the CO^ 
ice cap surface. 

These simulations do not directly require an accounting 
of incoming solar energy flux over all spectral wavelengths. 
The effect of this added energy source should have' the 
following repercussions: 

(1) a percentage of the total incoming solar flux will be 
conducted into the ice cap, causing corresponding tem- 
perature changes to occurj 

(2) a fraction will be’ used in the initiation of 

phase changes (solid C0 2 to gaseous C0 2 ) and the re- 
sultant C0 2 mass transfer by diffusion in the vertical 
will alter the vertical density distribution and also 
the thermal energetics; 

(3) a fractional amount of the incoming solar, radiation, 
will be attenuated through absorption and scattering, 
by the dust layer and may contribute to warming. 

The possibility remained that the Mariner 9 IRIS data was 
obtained over a frost - free region. To gain deeper insight 
into such a situation, another simulation was conducted 
with due consideration given to the interaction between 
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the ground and atmosphere. Pig. 13 depicts the dust- 
laden 1800 LMT initial temperature sounding adjusted 
to a surface temperature of 230K, together with the 



220 225 230 235 240 245 250 255 

Temperature, °K 


Fig. 13 Atmospheric and soil temperature profiles as 
a function of the hour (L^T) for a frost-free 
sand surface at 80° S; specifies the time 
of the initial sounding. 
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initial soil temperature profile extending to a depth 
of 60 cm. The results of a 24 hour simulation reveal a 
substantially-reduced overall layer cooling, except for 
the region above the dust layer. The magnitude of cooling 
above 8.5 km is markedly lower as can be seen by comparison 
with the soundings in Fig. 9. 

The effect of a thermally-conductive soil interacting 
with the lower polar atmosphere substantially alters the 
atmospheric thermal structure. A strong super-adiabatic 
layer extending from the surface to about 425 m is evident 
at 1130 LMT; by 1530 LMT this layer had expanded to the 
800 m level. The effect of the dust layer is seen by com- 
paring the 0400 LMT with the 1500 LMT soundings. Warming 
is observed below 7*5 km, while cooling occurs above. Be- 
tween 7*5 and 8.5 km, a super-adiabatic temperature dis- 
tribution is first noted at 1330 LMT . This layer remained 
super-adiabatic throughout the rest of the simulation due 
to the inability of the model to account, for convective 
adjustment. These simulative results support the fol- 
lowing tentative conclusions: 

(1) the condition of the lower boundary surface, 

soil versus solid CC> 2 , substantially alters the 
thermal energetics in the adjacent atmosphere. 
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(2) . this distinction, must be taken into consideration 
for any investigation of the vertical structure 
of temperature utilizing the integral inversion 
technique; for this effect may be as important 

to the derived inverted temperature profile as 
the presence of dust particulates, and 

(3) a more refined simulation model incorporating the 
effects of solid CO^ molecular thermal conduction, 
phase change, sublimation and vertical mass dif- 
fusion, should be developed for the polar frost 
cap region. 
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5.0 Theoretical Framework of the Generalized Problem 


5.1 Background & Orientation 

In the areas of non-LTE radiative transfer and local 
thermal energetics simulations, considerable developments 
have separately been achieved. This section of the report 
presents a detailed synthesis toward a rigorous remodeling 
program. 

The most direct impact of a dust storm on the radiation 
field of the atmospheric boundary layer is through 1) the 
existence of scattering and absorbing aerosol, 2) the close 
proximity of a lower boundary with specific radiative 
characteristics ,. and 3) the strong variation of the ab- 
sorptive and emissive powers at the top of the dust layer. 

The implication of 1) is that a rigorous modeling of radia- 
tive transfer is required which must include the formalism 
of non-conservative multiple Hie scattering by haze layers 
in the context of an inhomogeneous structure of atmospheric 
stratification. The interaction of the radiation field with 
the particulates as well as with the CO^ molecules controls 
the radiative heating/cooling rate. The mathematical treat- 
ment must consider polarization, in which case the equation 
of transfer is replaced by a matrix equation (Busbridge,1960) . 


54 



The solution of the matrix integro-differential equation 
of transfer has been discussed by e. g. , Adams & Kattawar 
(1970; "invariant imbedding" approach), Hanson (1 969 , 1971; 
"doubling method"), and Herman and Browning (1965; modified 
Gauss-Seidel interation). Bergstrom (1971) and Dave (1970) 
showed the iterative method to be sufficiently flexible. 

Shahrokhi and Wolf ( 1968 ) applied a finite difference 
scheme to the local integro-differential equation of radia- 
tive transfer. The Gauss-quadrature approximation of the 
source function integral yields a system of ordinary dif- 
ferential equations which lend themselves to finite dif- 
ference evaluation by direct iteration. 

The integration of the phase function and the scattering 
and extinction coefficients has most often been accomplished 
by using the modified gamma size distribution for the 
particulate polydispersion (Dave, 1969; Deirmenaj ian, 1969 ). 
Some information for fitting this distribution is available 
in the IRIS and other measurements of Mariner 9 . 

Construction of the atmospheric thermal energetics sim- 
ulation involves the calculation of radiative, conductive 
and convective heat flux divergences throughout the 
atmospheric layer structure, and the formulation of the 
balance of component fluxes at the interface between the 
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Martian soil and the overlying atmosphere. The resulting 
non-linear integro-dif ferential energy balance equation 
must then be solved for the temperature as a function of 
time and space. 


5.2 basic Definitions and Fundamental Relations of 
Radiative Transfer 


For detailed derivations of fundamental radiative transfer 
equations, we refer to the standard works (Chandrasekhar, I960 
Preisendorfer , 1965). In this section, the matrix relations 
will simply be listed, with some relevant remarks. 


Sekera ( 1963 ) has provided the local form of the equation 
of radiative transfer appropriate to the prototype problem: * 


iK V 


-[k 

+ K* U)X S u, u) 


( 26 ) 


with 


exp(--c^*) 

-*1 -UT 

+ w I (27) 




*Bold-face letters denote column matrix (e.g.,I ), and 
underscored letters denote 4x4 matrices (e.g., v jP)» 
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In the above: H :Matrix of Stokes parameters 

representing the diffuse radi- 
ation field for depths, zenith 
angle arc and azimuth <p \ 

^ :mass density 

:mass absorption coefficient 
:mass scattering coefficient 

(3 rmatrix of Stokes parameters re- 
presenting the unpolarized true 
emission ; 

{p rnormalized phase matrix for the 
direction pair (/*,<$>) m'j#’) • 

Diermendj ian ( 1 S 6 9 ) and Samuelson ( 1965 ) have shown that 
the phase matrix may be factored into a product of a scalar 
function depending on the size spectrum and number density 
of scattering centers, and a matrix containing only angular 
dependence. Similarly, the appropriate mass absorption and 
scattering coefficients may be derived for a heterogeneous 
polydispersion by a suitable averaging process. It is 
assumed that such operations have been performed in (26), 
(27), and subsequent relations. Introducing the total thick- 
ness, 

i 
00 

and albedo for single scattering. 


( 28 ) 


( 29 ) 


57 



(30) 


(26) and (27) may be recast as 

1 (Y " 

+ 4Tt J J S 'A 

“• a J * 

In contrast to the local form of the radiative transfer 
equation, global formulations have been developed resting 
on optical properties of a finite mass of atmospheric layer 
material. Conventionally, transmission and scattering mat- 
rices are defined (Chandrasekhar, I960; Sekera, 1963* 1966). 
The global point of view, which proved to be very successful 
in treating astrophysical systems, is only of limited value 
in the problem area of interface and boundary layer transfer. 
Therefore, its theoretical formulation and computational 
treatment will not be presented here. 

For the intensity matrix, boundary conditions must be de- 
fined which are given here in the standard form: 

= 1*1 d)-/k ,4>) - V 

-jr „ (3D 

Iv/CV,/,?) = TE v (t Ti >A,tf) = <D 

The asterisk denotes quantities referring to illumination 
from below. As explained by Chandrasekhar (i 960 ), solu- 
tions to the "planetary problem" (i.e., that with a re- 
flecting lower boundary) can be expressed in terms of the 
solution to the standard problem, by linear superposition 
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and slight modification of the_ matrices . Similarly, the 
effect of true emission at the lower boundary may be ac- 
comodated by adding an appropriate unpolarized isotropic 
intensity component to the upward intensity at the lower 
surface . 


5.3 Solution of the Local Radiative Transfer Equation 

An examination of the local radiative transfer equation 
reveals basically two solution approaches, depending on 
the choice of the dependent variable. If we choose the 
spectral specific intensity, then (26) is an integro- 
differential equation in which the angular variables appear 
as parameters. On the other hand, we can use the transfer 
equation to eliminate the intensity matrix from the source 
function definition ( 27 ) in which case the so-called aux- 
iliary equation is generated. In the latter relation, the 
source function plays the role of dependent variable. 

The auxiliary equation is of considerable theoretical 
utility, particularly for existence and uniqueness proofs, 
but has found computational use only in certain restricted 
contexts. These special cases are beset with many of the 

difficulties encpuntered with the global methods, because 
they were developed primarily for isotropic or Rayleigh 
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phase matrices. For these simplified scattering problems, 
the governing auxiliary equation becomes the familiar Milne 
equation (or, a simple variant of this form). The solution 
technique usually employed is the Neumann series (Busbridge, 
I960; Dave, 1964, 1965 ; Irvine, 1968 ; Preisendorfer, 1965 ). 

Many results have been obtained for multiple Rayleigh 
scattering problems, for which it can be shown that higher 
order terms in the Neumann series correspond to successively 
higher orders of scattering (Dave, 1964). As might be ex- 
pected, the extension of these analyses to the context of 
non-conservative Mie scattering is not at all apparent. 
Furthermore, convergence of the successive approximations 
is very slow for moderate optical thicknesses. 

Herman and Browning (1965) were among the first to regen- 
erate interest in the fundamental local radiative transfer 
equation, and appeal to a well-known iterative approach for 
its solution. It represents a true departure from the re- 
strictive confines of the astrophysical setting, is well- 
suited for simulation on a digital or analog computer, and 
yields directly the vertical distribution of spectral specific 
intensity, rather than the reflection and transmission 
matrices for a single optical depth. In principle, the ap- 
proach is suitable for an arbitrary phase matrix, and vert- 
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ical variation of scattering albedo. Thus, there are fewer 
a priori restrictive assumptions, such as the required sep- 
arability of the phase matrix or vertical homogeneity 
(e.g., Chandrasekhar, I960; Sekera, 1963, 1966). 


Recognizing the potential power of this basic approach, 
Dave (1970a, b) has developed an approximate analytical 
formalism for the Mle phase matrix, aimed at decreasing 
computing time with comparatively little sacrifice in 
solution accuracy. 


The local transfer equation in the form 

♦> ‘ ' 

'■? '•* 1 J/ V J -If. . un. )] ip it. ) 


(30) 


may be formally integrated with respect to optical depth, 
yielding a non-homogeneous matrix Fredholm integral equation 
for the Stokes parameters representing the spectral specific 
intensity. We shall, in what follows, suppress the true 
emission term in the matrix source function. Inclusion of 
this term creates virtually no additional difficulties for 
the solution algorithm, the Planck function simply repre- 
senting an additional term in the non-homogeneous portion of 
the integral equation. By the same token, we may assume 
that the single scattering albedo is unity, without any loss 
of generality. In this case, the equations governing the pth 
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Stokes parameter of the upward and downward intensities at 
the level T are 

+ c^ii) (32) 

Tt 

= [ \ -P >/,*') olio’ 

+ (33) 

In deriving the preceding relations, it has been assumed 
that the upward intensity at r«r T and the downward intensity 
at X-o are each identically zero. (Again, it should be 
apparent that there is no loss of generality). Thus, half 
of the intensities at r=r T and 'x-o are known, and half are un- 
known. 

Relations (32) and (33) constitute a system of simultaneous 
equations for the Stokes parameters representing the intensity 
at the level X . Kerman and Browning’s approach to the 
system can be described as follows. First, it is assumed 
that the dependent variables are known at some level Y n 
Then, the intensity at some distance from the level is 
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given by 1 (we restrict ourselv.es temporarily to the. upward 
intensity ) : ■ 




J /y 




-XlJyU. 


(34) 


where Tpl/^) is the value of the pth stokes parameter of 
the source function , -averaged over the increment hrt 
Writing the source function explicitly, and performing the 
integration indicated in (34), we obtain 


V A7f 


~f ^ = + O-e^l [ [f X 4 ‘ yne'd&J*' 

* ? n W) : 1^6)1,- ] 


(35) 


where X<^ and are the mean values of and over the in 
terval MT . Tnese averages are taken to be equal to the 
quantities evaluated at the midpoint of the layer. 


Finally, the angular integrations are approximated by sums 
over a discrete set of angular increments, and average values 
extracted. Performing all of these operations, and then 
solving' for. the intensity at the (n+l)st level, we get (Herman 
and Browning, 1965, eq.2.12): - . 


T w«0 A -z &C/14. n 






4 (l ‘ e 


(36) 
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The numerical procedure is initiated by writing the pre- 
ceding relation, for the first level below the "top" of 
the model atmosphere. To evaluate the summation term, 
values of for zenith angles corresponding to upward 
intensity are required. Since these are unknown, they 
are estimated for this first iteration to be zero. Vie 
may now use ( 36 ) to estimate T ^ (/* t <p ) for all p and n, 
down to the lower boundary, thus yielding the first iterate 
for the upward intensities. 

Now, a relation analogous to (36) may easily be derived 
for the downward intensities. The first iterates for these 
quantities are obtained by beginning at the lower boundary. 
However, this time, an estimate for the intensities required 
for the summation term is available from the first iteration 
for the upward intensities. In this manner, we obtain the 
first iterate for the downward intensities. The procedure 
from this point on is evident. The new estimates of upward 
intensities are used in (36), thus generating the second it- 
erate for upward intensity. The iteration is continued until 
successive values of intensity at each level agree to some 
preassigned tolerance. 


5.4 Solution of the Thermal Energetics Problem. 

For discussion purposes, the following pair of equations 
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may be taken as those governing thermal energy exchange 
in a dust-laden gaseous system (neglecting phase changes): 


( S.^)TT ic,5 4.1 = f s',-b) AlJ 

v ’ Hit JJ 

+ kj e ® v lT 


(37) 


? C T 




to 

+ 1 K V il «W- ■»» ® v lTir,t j]J }j„ 


(38) 


The first relation represents the non-LTE radiative transfer 
equation (see Sect. 5»2), and the second relates the local 
time rate-of-change ~of temperature to the divergence of 
the total heat flux. K(r ,g- # ,t) is a generalized heat dif- 
fusion coefficient which would include molecular-conductive 
and/or moderate turbulent heat exchange mechanisms as a 
function of position r and the static stability C*, The 
symbol £(.)][ designates the sum of the first 2 Stokes matrix 
components. The unit vector s_ indicates the direction in 
which the radiation field is considered. 

It should be noted here that it is the radiative net flux 
divergence which determines the radiative portion of heating. 
The second term on the right-hand side of eq. (38) represents 
a transformed expression for the radiative net flux divergence 
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as presented by, e.g., Bergstrom (1971) and Kondratyev 
(1969). 

For notational convenience, the relations (37) and (38) 
may be cast Into the form 


(39) 

and 

(1)0) 

where £ ( ) Is a linear integrodifferential operator, 

$ ( ) Is a linear parabolic operator, and H ( ), and G ( ) 

are linear and nonlinear integral operators, respectively, 
whose explicit form may be inferred from relations (37) and 
(38). 

Among the computational algorithms applicable to the non- 
linear functional relations (39) and (40) are finite differ- 
ences schemes (Ames, 1970), nonlinear optimization techniques 
(Beltrami, 1970), numerical integral transform methods 
(Bellman, 1969), and linearization approaches (Viskanta and 
Grosh, 1962; Lick, 1964). It appears, however, that the most 
promising approach will be through the use of the Newton- 
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Raphson-Kantorovic (N-R-K) iteration (Lee, 1968 ). 

The N-R-X-iteration is a function-analytic general- 
ization of the Newton root-finding .scheme for nonlinear 
scalar equations. It replaces the nonlinear functional 
equation with a sequence of linear relations whose solution 
converges to the solution of the original problem, under 
certain restrictions on the functional operators and initial 
iterate. The iteration scheme for relations (39) and (40) 
would in this context take the form 

jto - IWlb 

and 

(1|2) 

where ( ) f represents Frechet-dif ferentiation (Rail, 1969). 
The eqs . (41) and (42) constitute a linear system for deter- 
mination of the (n+ l)st iterates of temperature and radi- 
ative intensity, which may be solved by the methods employed 
by Lee ( 1 9 6 8 ) . 
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The Newton-Raphson-Kantorovic method generally ex- 
hibits quadratic convergence behavior, resulting in a 
more limited number of iterations required for a given 
accuracy, as compared with other methods, e.g., con- 
traction mapping. This advantage is particularly attractive 
when large gradients of the solution vector may be expected 
near the boundaries of the simulated system (Bellman and 
Kalaba, 1965 ). Such is the case for temperature near the 
surface or at the top of dust layers. 
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p .5 Modeling of the Time-dependent Thermal Convection 
Layer Induced by Radiative Surface Heating. 

While the scale analysis of Gierasch and Goody ( 1 96 8 ) and 
calculations of relaxation times for the Martian troposphere 
(Goody and Belton, 1967) indicate that advection of tempera- 
ture by tne planetary-scale wind field is of less importance 
in determining the vertical temperature profile than the 
effects of radiative transfer, it does not follow that the 
effects of small-scale convective heat transport are neg- 
ligible. The radiative equilibrium calculations of Prab- 
hakara and Hogan ( 1 965 ) > and time-dependent studies by 
Gierasch and Goody ( 1 968 ) and Pallmann and Dannevik ( 1972 ) 
indicate that stat ically-unstable thermal stratification is 
induced in the lower scale height by the intense surface 
radiative heating in regions away from the polar cap. Such 
instability would lead to development of a free-convect ion 
layer whose depth would vary over the daylight hours. In 
addition, it is possible that a detached convection layer 
would develop at the top of an optically-thick dust layer 
as a result of efficient radiative heating within the dust 
layer- combined with a sharp variation of emissive power 
across the top of the layer. 

An important feature of such thermal convection layers 
is that' they are capped above and/or below by regions which 
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are stably stratified, and thus depend upon penetrative 
mechanisms (convective heat flux opposite the mean temper- 
ature gradient) for their growth and maintenance. It is 
difficult to model this phenomenon with an Austauch 
parameterization of convective heat flux, because the 
latter requires a negative eddy conductivity to describe 
penetrative convection (Kuo, 1968). Furthermore, the 
application of similarity theory usually requires an as- 
sumption relating to existence of a constant flux layer. 

An alternative approach to numerical modeling of the time- 
dependent free-con vection layer is to calculate explicitely 
the convective heat flux by determining the dynamical 
motion and temperature fluctuation fields. This requires 
solution of the equations of motion, similar in some 
respects to that undertaken in the nonlinear Rayleig'n- 
Benard problem. However, for insolation to provide a 
natural drive for the convection layer, a thermal lower 
boundary condition is required which is considerably more 
general than that which is conventionally applied to the 
classical Benara problem. The required thermal boundary 
condition is, in fact, the same which has been applied to 
the radiative-conductive model described in this report. 

If the maximum depth of the surface-based free-con vection 
layer is small compared with the scale height of the isen- 
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tropic atmosphere and we are willing to exclude acoustic 
waves from the analysis, the appropriate equations for the 
required numerical model are the Boussinesq approximations 
to the Navier-Stokes , energy, and continuity equations. 
Herring (1963, 1964) has developed a spectral model based 
on these relations, employing the "mean-field approximation. 
In this formulation, it is assumed that the nonlinear inter- 
actions among fluctuations of temperature and velocity are 
of lesser importance in determining the vertical convective 
heat transport than the interactions between vertical velo- 
city fluctuations and the mean (horizontally-averaged) vert- 
ical temperature profile. Thus, the model is capable of 
simulating the important mechanism of penetrative convection 
but is considerably - simpler to handle numerically than the 
system Including the fluctuating interactions. 

Musman (1968) has shown that Herrings’ model is capable 
of reproducing in some detail the circulation induced by 
a steady penetrative convection and it also appears 
(Elder, 1969) that the model is quite successful in the time 
dependent case. 

In the mean-field Boussinesq approximation, the system 
of governing equations is closed without the necessity of 
computing the horizontal velocity components and the de- 
viation of pressure from hydrostatic equilibrium. Thus, 
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there are three equations for determination of the vertical 
velocity fluctuations, horizontally-averaged vertical temper- 
ature profile, and the fluctuation of the temperature field 
from the latter profile. In the spectral version of this 
model, the boundary conditions are automatically satisfied 
because of the special nature of the basis functions used in 
the series expansion. 

Herring’s model can be extended and modified to be com- 
patible with the radiative-conductive model outlined in 
this report. The required generalized thermal boundary con- 
dition would still be applied to the subsoil domain, to 
determine the new surface temperature at the end of a time- 
step. This temperature would then furnish the lower boundary 
condition for the spectral convection model, which would 
include the effects of radiative transfer. This model would 
then generate a new atmospheric temperature profile, based 
on radiative, conductive and convective heat transport. The 
new profile would furnish the downwardly-directed radiative 
fluxes required for calculation of a new soil surface temper- 
ature and the cycle would be repeated. In this manner, the 
time-history of the subsoil and atmospheric temperature pro- 
files can be computed, including the combined effects of 
radiative, molecular-conductive, and thermal-convective heat 
transport mechanisms. 
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.6.0 Summary of Research and Major Findings 

As described in Sections 1 and 2, the wealth of in- 
formation provided by the Mariner spacecraft series has 
generated an adequate basis to allow inferences into the 
primary physical mechanisms likely to be active in Martian 
atmospheric processes. The high concentration of radi- 
atively-active constituent and the relative tenuity of 
the Martian atmosphere result in an intense radiative 
response capability, such that tropospheric heat exchange 
at scales greater than 10 m is dominated by radiative 
transfer. By contrast , terrestrial radiative relaxation 
of thermal perturbations is at least 30 times slower, 
i.e., ^ 10 days for large-scale motion. In addition, 
the global extent of sand-like dry soil material 
permits insolation to induce a large diurnal soil surface 
temperature variation in the regions away from the polar 
ice cap. Over the dry ice, the surface temperature is 
fixed at the C0 2 sublimation point, regardless of the 
variation of insolation. Because of its radiative re- 
sponsiveness, the lowest several km of the Martian at- 
mosphere directly responds, in either case, to the thermal 
conditions at the lower boundary. The resulting large 
thermal gradients near the ground-atmosphere interface 
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make molecular thermal conduction an important mechanism 
in coupling the sub-surface and atmospheric thermal struc- 
ture. When atmospheric strata are dust-laden, an ad- 
ditional radiation table is provided at the top of the 
dust layer, which acts radiatively like a "secondary sur- 
face," significantly influencing atmospheric thermal struc- 
ture in the layers above and below. Thus, secondary thermal 
boundary layers are induced. 

The physical mechanisms having primary roles in the deter- 
mination of thermal structure in the Martian ground-atmo- 
sphere system are such that a numerical model of the time- 
dependent structure of thermal boundary layers can be con- 
structed without the need to consider directly the atmos- 
pheric motion field. This model has been described in 
Sections 2 and 3. Some of the more important results de- 
rived through basic numerical experiments with the model 
are as follows: 

(1) when there is little dust present in , the Martian 
atmosphere, and in regions away from the polar 
ice caps, weak-and strong-line radiative com- 
munication produces three primary atmospheric 
regimes. The lowest, characteristically of 3 km 
depth, responds directly through strong-line 
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radiative transfer and molecular thermal 
conduction to surface temperature variation. 

The middle regime, extending to~30 km altitude, 
is in weak-line communication with surface temper- 
ature variation, but also undergoes direct solar 
heating and losses to outer space. The upper 
regime is dominated by the upper boundary con- 
ditions, i.e., solar heating and emissive losses, 
and is affected little by conditions at the sur- 
face . 

(2) the presence of an atmospheric dust layer which 
is optically-thick in the near-IR ranges induces 
secondary thermal boundary layers at the ex- 
tremities of the layer, because of the sharp 
variation of radiative properties. For example, 
diurnal temperature variation at the top of a 
dust layer of 30 km thickness may be as much as 
50 $ of the soil-surface temperature variation. 

(3) midlat itudinal temperature soundings derived by 
inversion of Mariner 9 IRIS measurements obtained 
during the dust-laden period of the mission, under- 
go considerable modification when used as initial 
data in the simulation model. There is some in- 
dication that the IRIS-derived profiles, obtained 
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from an inversion technique assuming a dust- 
free atmosphere, represent a highly smoothed 
version of the actual thermal structure of the 
dust-laden layers. 

(4) in the polar region, the temperature of the dry 
ice surface is rigidly maintained at 148K, the 
sublimation temperature of CO,, at a Martian sur- 
face pressure of 6 mb. The heat transfer mech- 
anisms at such an interface are dominated by this 
constraint. The simulative outputs indicate that 
the planetary atmospheric counter-radiation plays 
an important role in the heat flux balance at the 
interface. The persistency of the heat sink at 
the dry ice surface removes rather large amounts 
of heat from the atmosphere by conductive and 
radiative processes. In contrast, polar ice-free 
sand surfaces are not only warmer and maintain 
variable temperatures but also permit the tropo- 
spheric , layers to respond to surface heating/cooling. 

(5) the time-dependent simulated temperature structure 
in the atmosphere over a frost-free surface of the 
polar region is found to vary between 220K and 
235h. - These results are rather closely represented 
by the Mariner 9 IRIS sounding. However, if the 
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underlying surface is frost-covered the simulated 

0 

temperature profiles are substantially cooler. 

The outputs from our model indicate that a large 
fraction of the frost-surface emission does not 
penetrate the dust layer. 

In Section 5» the theoretical framework of a refined 
version of the simulation model has been presented. Ex- 
perience with the current version of the model shows that 
among the areas requiring refinement are (1) the treatment 
of the impact of non-conservative multiple Mie-scattering 
by particulates on thermal energetics transfer, and (2) the 
incorporation of free-convective heat transport mechanisms 
to accomodate the impact of st atically-unstable temperature 
distributions induced by radiative exchange. 

It is shown that the impact of multiple Mie scattering 
on thermal energetics can be incorporated by solving sim- 
ultaneously the heating-rate equation and the radiative 
transfer equation including the scattering terms. A local 
iterative scheme based on the Newton-Raphson-Kantorovic 
and Gauss-Seidel algorithms is suitable for the resulting 
system of equations. 
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In addition, a method for incorporating free-convectiwe 
heat transfer into the simulation model has been out- 
lined. This scheme is based on a spectral version of the 
equations applicable to the Rayleigh-Benard problem, with 
a generalized lower boundary condition to allow thermal 
coupling between the ground and atmospheric subsystems, 
and to permit insolation to provide a natural driving 
mechanism. It is shown that this approach is compatible 
with the radiative-conductive model including the effects 
of particulate multiple scattering and absorption. 
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